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acceleration due to gravity; m/s 2 
height of liquid (vapor); m 
gas/vapor or liquid mass flow rate; kg/s 
mass of a control volume; kg 
(partial) pressure; Pa 

heat flow rate from or to a control volume; W 
heat flow flux from or to a control volume; W/m 2 
absolute temperature; K 
time; s 
volume; m 3 

velocity of liquid or vapor; m/s 
power on or by a control volume; W 
density; kg/m 3 

constant pressure (volume) specific heat; J/K-kg 

dimensionless Grashof, Nusselt, Rayleigh, Reynolds, and Prandtl numbers 
specific enthalpy; J/kg 

hydrogen and helium gas constants; J/(K-kg) 

specific internal energy; J/kg 

convection heat transfer coefficient; W/(m 2 -K) 

coefficient of volumetric expansion; 1/K 

ratio of specific heats; cp/cy 

thickness of boundary layer; m 

thermal conductivity; W/(m-K) 

dynamic viscosity; kg/m-s 

kinematic viscosity; m 2 /s 
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A = lateral surface of boundary layer element; m 2 
/ = dimensionless resistance coefficient 
H = height of external tank; m 

R = radius of external tank, storage tank, or a pipe; m 
S = cross-sectional area of vent valve area; m 2 

x = vertical position counted with respect to bottom of corresponding control vlume; m 

A = thickness of the horizontal layer; m 

A = dimensionless vent valve position 

Subscripts 

a = ambient 

B = bulk layer element 

C = critical point 

e = external heat or mass flows across the control volume boundaries other than liquid/vapor interface 
/ = gaseous hydrogen vapor film 

g = gaseous helium 
hor = a horizontal surface 

k = total number of liquid control volume layers 

L = boundary layer element 

l = liquid hydrogen 

lam. = laminar flow 

Is = saturated liquid 

Iv = liquid/vapor phase transition 

n = total number of layers in vapor control volume 

S = liquid/vapor interface 

turb = turbulent flow 

v = gaseous hydrogen vapor 

vs = saturated vapor 

w = tank wall 
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Superscripts 

cond = conduction 

conv = convection 

A = gaseous hydrogen or helium 


I. Introduction 

Recently, a dynamical model was developed for a complex spatially-distributed system of liquid hydro- 
gen (LH2) loading that involves the storage (ST) and external (ET) tanks as well as the transfer line [1], This 
reduced model is based on a set of coupled ordinary integro-differential equations, which is shown to be 
well-suited to describing a generic cryogenic loading system, including one used for Space Shuttle fueling 
in both the nominal and major faulty regimes. The model accounts for the pressurizing helium gas injected 
into the ET and can be easily modified for similar cryogenic propellant loading systems. 

Analysis has been performed for both the nominal regime and the effects of several primary faults, such 
as gas leaks in the ullage space of both the ST and ET, as well as clogging of the ET vent valve on the history 
of both tanks [1,2], It is found that each of the faults is characterized by quite a pronounced dynamics. The 
most interesting observation is that some of the faults, such as a substantial clogging of the ET vent valve, 
usually have a slight effect on the integral variables, namely on the volumes of LH2 in the tanks or even on 
the ullage temperatures and pressures. Yet, the dynamics of the ET vent valve switching that causes pressure 
oscillations are shown to be extremely sensitive even to small deviations from the nominal regime and can 
be used for early fault identification by means of real time sensor data analysis. (The valve in question works 
between two threshold positions in order to maintain a required ullage pressure [1].) It means that one can 
not only detect, isolate and identify different faults based on the analysis of the sensor data characterizing 
the filling dynamics, but also infer the parameters of the model to be further used for prediction of future 
behavior of the system [2], 

However, the above model, allowing for a fairly adequate description of the ET ullage history [1,2], does 
not account for a substantial temperature stratification that may happen [3-6] in both the ullage [7-10] and the 
liquid [9, 10] control volumes of the external and storage cryogenic tanks. The major goal of the present paper 
is to develop a model that correctly describes possible thermal stratification in both the ullage and liquid bulk 
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parts of a tank that is based on ordinary differential equations, similar to those discussed in [1], short of a full- 
scale computational fluid dynamics (CFD) modeling scheme [3-8], so that it can be implemented efficiently 
and used for online state estimation, fault diagnosis, and fault prognosis. A distinctive feature of this model 
is that it takes into consideration non-equilibrium mass and energy exchange processes of condensation- 
evaporation at the liquid-vapor interface between the liquid and gas, and the thermal convection in both the 
liquid and vapor parts of the tank, which is driven by complex phenomena developing in the boundary layers 
adjacent to the tank walls as well as to the liquid-vapor interface. 

The temperature stratification effects under consideration may be responsible for the possible appearance 
of the overheated layers [11, 12] of cryogenic fuel in the external tanks, and would be crucial to designing 
effective pressure and temperature control in the storage tanks [7-10, 13]. The model is meant to be ap- 
plied towards treating very important aspects of cryogenic fuel management technology [13], such as the 
minimization of the boil-off effects and the optimization of fuel bias, i.e. of the excess of fuel load above 
nominal aimed at preventing the cryogenic fuel from overheating due to radiation from the hot jet and heat 
conduction transfer through the ET walls from the engine, thrust cone, and common bulkhead. The over- 
heated cryogenic liquid on its way to the engine may reach a pump impeller thus leading to the development 
of violent cavitation [14] in the areas where pressure is less than the saturated pressure of the overheated fuel. 

The contributions of this work are as follows. We develop a general model of a cryogenic tank including 
temperature stratification effects in the liquid, vapor, and tank control volumes, for an arbitrary number of 
lumped-parameter layers. We develop a simulation that implements the model and allows the investigation 
of the dynamical behavior under both storage, loading, and unloading regimes. In particular, we show that 
the presence of a vent valve in the tank used for the purposes of maintaining pressure creates temperature 
waves that propagate through the layers. Simulation results are presented and detailed analysis is provided. 

The paper is organized as follows. Section II describes the cryogenic tank model, which incorporates 
temperature stratification effects driven by convection in both the liquid and vapor parts of the tank. Sec- 
tion III deals with major results of our simulations of the LH2 cryogenic tank in different performance 
regimes. Section IV concludes the paper. 
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(a) Side view. (b) Top view. 

Fig. 1 Mass flow in a propellant tank with three control volumes (vapor, film, and liquid). 

II. Description of the Model 

Consider a cylindrical tank of radius R where the gravity/inertia force is parallel to its axis. It is partially 
filled with liquid hydrogen (LH2); its ullage volume is filled by a mixture of gaseous hydrogen (GH2) and 
gaseous helium (GHe). The liquid and ullage (vapor) control volumes (CVs) are separated by a thin massless 
film of saturated hydrogen [1] (see Fig. 1). The following subsections describe the liquid and vapor CVs and 
the heat and mass flows depicted in Fig. 1. 

A. Vapor Control Volume 

Divide the vapor (ullage) CV into n > 3 layers perpendicular to the axis of the tank. The lowest one 
( i = 1) is adjacent to the saturated film CV and exchanges mass and energy with it. The upper one (i - n) 
is at the top of the tank and exchanges heat with the top wall of the tank, as well as mass and energy with the 
environment through the vent valve, as shown in Fig. 1. Any internal horizontal slice (i = 2, . . . , n — 1), in 
turn, consists of two sub volumes: (i) a torus-shaped element of the boundary layer (L), which is adjacent to 
the lateral wall of the tank and exchanges heat with it; and (if) the cylinder-shaped bulk element (B), which 
is coaxial with the tank vertical axis x. Thus the total number of the vapor CVs is equal to 2 (n — 2) + 2 = 
2(n — 1). 

In the presence of the two gas components in the vapor CV: GH2 ( v ) and GHE (g), for the zth internal 


6 



vapor bulk ( B ) element (see Fig. 1), the mass and the energy conservation, respectively, yield [15-21] 


^(A) _ 7 (A) _ 7 (A) _ t(A) 

" L i,B ~ J i,B J i,BL 


( 1 ) 


- V V J W h (X) - ( J {x) + J (x) ) h (x) 

^ ~ vv i,B ^ 2-^ J i+l,B ri i+l,B \ J i,B ^ J i,BLJ fl i,B 

X —v,g A =v,g 


( 2 ) 


where A = v, g. Here, when the tank walls are hotter than the gas/vapor (T wv > T'v)i b), the flow vertically 
enters the element (i. B ) from the above layer (i + 1, B) and exits to an underlying bulk element (i — 1, B)\ 
the remaining part (i, BL) spreads radially outwards and enters the /th internal boundary-layer ( L ) element 
(see Fig. 1). 

Similarly, for the zth internal boundary -layer (L) element (see Fig. 1) 


• (A) _ hA) _ t(A) 7 (A) 
" l i,L ~~ J i~l,L J i,L ' J i,BL 


( 3 ) 


- V m (X) u [X) - 0 (v) - W {v) + V J (A) h (X) - J (X) h W + J (A) h () 

1. nl i,L a i,L — b!i,e yv i,L ' J i-l,L n i-l,L J i,L u i,L ' J i,BL n i, 


(A) ; (A) 


/(A) t (A) . t-(A) 7 (A) 


dt 

X—v,g A =v,g 

Here, the specific enthalpies and energies of the elements are proportional to their absolute temperatures: 


( 4 ) 


7 V(g) rji 

tl i,L{B) ~~ C P,v(g) 1 i,L(B) 


v{g) rji 

U i,L(B) ~ C V,v(g') 1 i,L(B) 


( 5 ) 

(6) 


For the zth layer, the power associated with the change in its volume, Vi^B(L) [1, 15]: 


T A v ) 


W, 


,B(L) ~ (Pv + Pa)^i,B(L) 


( 7 ) 


The rate of the heat transfer from the wall to the vapor boundary layer element [15-21] 


Q^=a z A z (T w>v -Ti$) 


( 8 ) 


where A; t = 27rf?A, is the lateral surface of the boundary-layer element of the height A,;. (In this paper the 
tank wall in contact with vapor is treated as a separate CV with temperature T wv compared with the tank part 
in contact with LH2 which has temperature T w i.) 

In the case of the isothermal vertical wall (Fig. 2), the heat transfer coefficient can be calculated as [7, 8, 


16-21] 


K 

OCi — 

Xi 


0.68 + 0.503(7?^ T-) 1 / 4 , 
0.15(.Ratf) 1/3 , 


10 5 < Ra < 10 9 
10 9 <Ra< 10 11 


( 9 ) 
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and the thickness of the velocity boundary layer, within which the vertical component of the velocity is not 
very small (Fig. 2), can be calculated as [16-21] 


3.93 


Si — Xi \ 


0.952 + Pr\ 1/4 
Gr;Pr 2 ) '' 


lam 


1/10 


0.565 (i±^) /p r 8 / 15 ) turb 


(15) 


It can be seen that the vertical mass flow rate is fully defined by the parameters of the corresponding boundary 
layer element (i, L). In the case of natural convection, the mass flow is treated as laminar when the Raleigh 
number (11), which here plays a role of the Re number, is small [15], usually Ra < 10 9 [15-20], The 
thickness of the thermal boundary layer, where the temperature changes from that of the wall T w to that of 


the bulk Ti t L 


Si,r = 8i/ Pr 1/2 


(16) 


is usually smaller than that of the velocity layer (Pr > 1) [16-20] (see Fig. 2). In this paper the effects of film 
boiling are ignored [14]. A slight modification of the above formulas, describing the boundary layer effects, 
is required if the wall is not considered isothermal (see [15-21] and Appendix A). 

For the lowest horizontal vapor layer (i = 1), which is in contact with the film CV, 


Av) _ T (v) 


"1,B 


= A v k - A V 1 - Jiv 


' 2,B °l,L 


Ad) _ j(,g) 


n,B 


( 9 ) 


= A , L 


d 

dt 


E = Qv- w^ B + e - jAIAA) - J ^ h v 


(17) 

(18) 
(19) 


A =v,g X=v,g 

Here, it is assumed that the lowest horizontal layer does not have the boundary element, so that the horizontal 
flow is equal to the upward vertical flow (1, L) and is given by (13). 

In accordance with [1] the rate of the heat exchange between the lowest horizontal gaseous layer and the 
vapor interface film CV, which is driven either by convection or by conduction, can be calculated as 


Q v =±m^{Q c v ond ,QT v } (20) 

('ll) 

where + is chosen when the temperature of the vapor-liquid interface, Tf = Ts < Tf B , and — otherwise; 
the rate of the heat exchange due to convection 


QT V = ^R 2 oAor ( T s - rg) H ( Ts - rg) 


( 21 ) 
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where the Heaviside function H takes care of the proper temperature condition required for convection; 
the heat transfer coefficient from a horizontal surface [16-20] 

K 0.54 f?c4 /4 , 10 4 < Ra R < 10 7 

Othor = < ( 22 ) 

I 0.15Ra]( 3 , 10 7 < Rclr < 10 11 

Here the Rclr number is given by (1 1) and (12) where T w — T t ,n should be substituted with T s — 1\js and 
x with R\ and 'k is given by (10). In (20) the rate of the heat exchange due to conduction, Q^ ond , when it 
dominates heat flow through the interface, is defined by the time evolution of the temperature of the interface 
(film) Ts(t) . All the details of computation of Q% md can be found in our previous paper [1], 

For the upper horizontal layer, which is in contact with the top of the tank and presumably with the vent 
valve, 


t tW I 7 

m n,B ~ Jn-1 ,L J n ,B + 


(23) 


J t E m i% U ni=Qto p -Wi%+ E (24) 

A =v,g X=v,g 

where the rate of the heat exchange with the top of the tank is given by 


Qtop = nR 2 a^ r (t wv - T™) (25) 

(v) 

with the heat transfer coefficient given by (22) where Ts — T) b should be substituted with T wv — B . 
Here, it is assumed that the upper horizontal layer does not have the boundary element (see Fig. 1). 


B. Liquid Control Volume 

Divide the liquid CV into k horizontal layers, similar to the vapor CV. (The total number of the liquid 
CVs is equal to 2 (k — 1).) Then for the jth bulk element (j = 2, . . , , k — 1); the levels are counted from the 
bottom of the tank up to the liquid/vapor interface) 


d 

dt 


r (,) 


(l) ( l ) 

m j,B U j,B 


__ t (0 _ t ( 0 j{l) 

~ J j + 1,B J j,B J j,BL 


— jF) hd) 

— vy j,B + J j + l,B n j + l,B 


J- l i + J ) 




(0 


j,BL 


h® 

n j,B 


(26) 

(27) 
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For the yth internal boundary-layer (L) liquid element, the mass and energy conservation, respectively, yield 
(see Fig. 1) 


fh^L = 


3 - 1 , l 


- 4'i + J 


( i ) 

j,BL 


dt m j,L U j,L - Wj,e VV j.L + J j-l,L n j-l,L J j,L n j,L + J j,BL n j,B 


Here the specific enthalpy and energy are close to each other and both approximately equal to [1]: 


(28) 

(29) 


h j,L(B) ~ u j,L(B) « c I T j,L(B ) 


(30) 


For the upper horizontal layer, which is in contact with the film CV (interface). 


• (0 t(0 I t(0 I T 

1 TO (0 u (0 —n _ j(0 ^(0 I j(0 ^(0 I j u 

d -j- " l k,B a k,b — vv k,B J k,B u k,b ^ J k-l,L n k-l,L ' J lv'Hs 


(31) 

(32) 


Here it is assumed that the lowest horizontal layer does not have the boundary element, so that the upward 
vertical flow is equal to the horizontal flow and is given by (13). 

Here, the rate of the heat exchange with the film CV 


Qi =±rnax{Qr d ,0r m ’} (33) 

where + is chosen when Tg > T^l , and — otherwise; the rate of the heat exchange due to convection 

Qi° nv = ( t kb - Ts) H (ig, - T S ) (34) 

with the heat transfer coefficient given by Eq. (22) with parameters corresponding to LH2 [15-21], 

For the lowest horizontal layer, which is in contact with the bottom of the tank, and, assumingly with a 
transfer line valve 


. (0 _ 


m ; = - J[ 1) B + J£> R - J L , 


t( 0 

J 2 


dt 


= Qbott - w<% - J[%h^ g + 4%h^ L - J lte (h^ B + V&2) 


where the rate of the heat exchange with the bottom 


(35) 

(36) 


Qbott = 7rR 2 a^ or (T wl - T^l) (37) 

with the heat transfer coefficient given by (22) where T wv — should be substituted with T w i — T^l- 
Here, it is assumed that the lowest liquid horizontal layer does not have the boundary element. (All the details 
of computation of the heat flow due to conduction, Q < f ond can be found in our previous paper [1].) 
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C. Summary of the Model 


Thus, dividing the liquid CV into k horizontal layers and the vapor CV into n layers, one 

deals with the following 2(5 n + 3fc) — 13 state variables: 4 (n — 1) + 2 (k — 1) = 2(2n + 

k - 3) masses m\' B , my' B , m\ ^ {L) , m\' B , [i = 2 - l;j = 2 ,...,k - 

1;A = v,g); vapor CV partial pressures and volume p v ,p g and V v ; 2 (n + k) — 1 temperatures 

Tf, T wU T wv ,T[% (i = 2, . . . , n - 1; j = 2, . . . t k - 1); and 4 (n - 2) + 

2 (k - 2) + 3 = 2(2 n + k) - 9 mass flowrates J j%L> J i% (* = 2 > ■ • • > n - 1; j = 

2, . . . , k - 1; A = v,g). 

There are 4 n + 2k - 3 constraints: 4 (n — 1) GH2 and GHE equations of state; 




-,(v) 


(38) 


where A = v,g; the film equation of state 


p v =p f (T f )=p c (T f /T c y 


(39) 


with pc = 1.315 MPa and Tc = 33.2 K and n = 5 being, respectively, critical pressure and temperature for 
hydrogen [1] (in (38), the partial pressures are considered uniform across the vapor CV, with the LH2 partial 
pressure equal to that of the saturated vapor film (39); 2 (k — 1) — 1 equations relating the liquid masses and 
volumes (because one liquid mass is fully defined by the total liquid mass, which, in turn is defined by (41)) 


in 


(i) 

iMB) 


P^Mbv 


(40) 


and an equation that relates the total tank volume to the GH2 volume and LH2 volume 


V = V v + Vi= | V™ + v^’l + £ (v$ + v t 


( V V 


(41) 


\ 7= 2 / \ j = 2 ) 

There are 6n+4fc— 4 ordinary differential rate equations for 2(5n+3/c) — 13— (4n+2fc— 9) = 6?r+4/c— 4 
independent variables: 4 (n — 1) GH2 and GHe mass conservation equations (1), (3), (17), (18), and (23); 
4(n— 1) GH2 and GHe energy conservation equations (2), (4), (19), and (24); 2(k— 1) LH2 mass conservation 
equations (26), (28), (31), and (35); 2 (k — 1) LH2 energy conservation equations (27), (29), (32), and (36); 
and two energy conservation equations for the two control volumes for the tank wall in contact with LH2 and 
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Fig. 3 Heat flow in a propellant tank with three control volumes (vapor, film, and liquid). 


the ullage, respectively (Fig. 3) 


C'wdirriyjiTwi ) / dt — 


Qwl Qe + C'vjdlwid'wl’) dl w i < 0 

Qwl Q e CwrilyjlTqjjv 1 dlyjl ^ 0 


(42) 


c w d(m wv T wv ) / dt — 


Qwv Qe + c w rn wv T w i i diyji < 0 


Qwv Qe + C w m VJV^UW 1 rilwl ^ 0 

Here the last terms in the right sides of Eqs. 42 and 43 describe the flow of internal energy along the tank 
walls due to motion of the liquid level and of the subsequent plane separating the two tank CVs (note that 
rhwi = — riiwv )• Also, in the above equations we ignore the direct heat exchange (due to conduction) between 
the two parts of the tank because the corresponding thermal length even for t - 10 3 s is shown to be much 
smaller than the size of those parts. The total heat transfer from the tank wall to the vapor CV (see (8) and 
(25)) 


(43) 


Qi v} = Q top + £ QS 

»= 2 

and the total heat transfer from the tank wall to the liquid CV 

k - 2 

Qe^ = Qbot.t + £ 
j=2 


(44) 


(45) 
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where 


Q g \ e = a jAj(T w i — T^g) (46) 

and Qbott is given by (37). Here the heat transfer coefficients a 3 are given by (8) with the parameters of 
liquid and A 3 = 2ttRAj is the lateral surface of the boundary-layer element of the height A j . 

And, finally, the condensation-evaporation mass flow rate through the LH2/GH2 interface, is related 
to the heat flows through the same interface, Q v (20) and Qi (33) by means of the energy conservation 
equation for the film CV. Ignoring the film mass, it can be written as 

Qv Ql T J lv ( h'vs — 0 (47) 

where the specific enthalpies are referred to those of the saturated vapor and liquid. The temperatures of the 
vapor liquid interface, Ts, is considered to be equal of that of the film Tf [1], 

III. Results 

The above described model of temperature stratification has been explored in an LH2 tank. The pa- 
rameter values of the cylindrical tank were chosen similar to that of the Space Shuttle ET [1, 2], provided 
in Appendix C. In Eqs. 42 and 43, the material parameters of the tank are chosen to that typical for alu- 
minum (see Table 1 in Appendix C). The tank is considered to be initially (at t = 0) half-filled with LH2: 
VJ(0) = V/2 with GH2 and GHe partial pressures p„(0) = p g ( 0) = 1 atm and initial temperatures of 
GH2/GHe mixture, LH2, and the tank, respectively equal to T v (0) = 71(0) = 20 K and T w (0) = 21 K. 
In addition, the tank absorbs heat from the environment (see Eqs. 42 and 43) at a rate of Q w = 20 kW in 
the blowdown regime. Throughout modeling, the number of horizontal layers in the liquid CV was set equal 
to that in the vapor CV, where n = k = 10, which yields 146 state variables and 95 ordinary differential 
equations. It is found that starting from n = k = 10 and on, such “integral” characteristics as the total 
ullage pressure p v + p g and the tank wall temperatures T w i and T wv converge and change no more than, 
respectively, 0.1% and 0.3% with the increase of the number of the horizontal layers. Those numbers can be 
taken as a measure of the uncertainty of the predicted results. 
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A. Blowdown Regime 


In the blowdown regime it is assumed that the ullage vent valve is closed ( J u ,e = Jg,e = 0), and LH2 is 
being removed from the tank at a constant rate of e = 10 kg/s. The case of normal gravity is considered 
first. Fig. 4 shows the masses and temperatures. (For this and all the subsequent figures higher layers are 
shown with lighter-colored lines.) LH2 mass (Fig. 4a) drops, due to the positive external mass flow, but also 
due to evaporation, while the GH2 mass (Fig. 4b) increases due to evaporation. The very same evaporation 
is responsible for cooling off the LH2/GH2 interface, and the liquid-vapor interface temperature Tg = Tf 
drops due to the evaporation ((Figs. 4c and 4e). The lumped-parameter ullage temperature is also shown in 
Fig. 4e as T v to demonstrate that the stratified temperatures behavior is consistent with the lumped-parameter 
formulation (see Appendix B). The expansion of the ullage volume leads to its initial cooling, which gives 
way to the subsequent ullage heating due to effective convection-assisted heat transfer from the tank wall 
(Fig. 4e). The temperatures of the liquid layers also increase due to heat transfer from the tank wall, as 
shown in Fig. 4c. Fig. 4d shows only the temperatures of the liquid bulk temperatures, since the stratification 
is not easily observed in Fig. 4c. Fig. 5 uses a heat map to show a visualization of temperature stratification 
in this case at different time instants. 

Fig. 6 demonstrates temperature evolution of different layers of the tank as it would be measured by 
sensors at fixed heights in the tank. Here, n = 100 is used to minimize an artificial “ripple” effect caused 
by transition from one layer to another that happens at a fixed sensor height. It can be seen that as the liquid 
level crosses the sensor at 12 m, the temperature experiences a sudden drop from the LH2 temperature to that 
of the GH2/GHe mixture. A small peak seen at the ullage temperature curves can be attributed to transient 
effects that develop at the early stages of the system evolution. 

Fig. 7 shows the masses and temperatures for the blowdown regime in elevated gravity (7 g). Here, 
the evaporation rate is somewhat larger than the case of normal gravity (compare Figs. 4b and 7b), and 
stratification seems to be reduced slightly due to the enhanced convection-driven heat exchange between 
the tank body and its contents (see Figs. 7c, 7d, and 7e). Fig. 8 shows the visualization of temperature 
stratification in this case. 

Fig. 9 shows behavior of the masses and temperatures for the tank at the conditions of microgravity 
(10 6 <7). Here, the heat exchange between the tank wall and its contents is suppressed, thus leading to 
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Fig. 4 Tank in a blowdown regime at normal gravity. 
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Fig. 5 Tank in a blowdown regime at normal gravity: temperature stratification. 


continuous cooling of the ullage due to its expansion. In this case, the suppression of the heat exchange, 
in contrast to the initial drop in the tank’s temperature T w i and T wv at normal (Fig. 4) and elevated gravity 
(Fig. 7), causes continuous heating of the tank wall and diminishes temperature stratification effects in both 
the liquid and vapor CVs. This is clearly shown visualization of the temperature stratification in Fig. 10. 

It should be mentioned here that in all cases, as expected, the temperature stratification is much more 
pronounced in the ullage CV compared to the liquid CV due to much more effective convection flow and 
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Fig. 6 Temperatures at fixed heights in a tank in a blowdown regime at normal gravity. 
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Fig. 7 Tank in a blowdown regime at elevated gravity 7 g. 


thermal conductivity in the LH2 (.//, = 5 kg/s in comparison with GH2/GHe where .7/ = 0.9 kg/s). Figs. 4e, 
7e, and 9e also demonstrate that the results of the simulations, which take into consideration the effects of 
temperature stratification, are fully consistent with previous ones (seen on T v curves) that are obtained in the 
framework of the lumped parameter method (see [1, 2] and Appendix B). 


B. Storage Regime 

The storage regime is now considered, where liquid is not added or removed from the tank ( e = 0). 
The tank, with initial temperature set at T w ( 0) = 20.3 K is absorbing heat at a rate of Q w = 10 kW that, for 


17 



Fig. 8 Tank in a blowdown regime at elevated gravity: temperature stratification. 
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Fig. 9 Tank in a blowdown regime at microgravity. 


the tank of a chosen size, roughly corresponds to the radiation heat flux from the Sun. Fig. 1 1 demonstrates 
temperature evolutions in the unvented tank (J v ^ e = J, he = 0) at the conditions of normal gravity. Here, 
evaporation is greatly reduced (Fig. 11a). Substantial temperature stratification in the ullage volume can be 
observed. Also, it can be seen that in contrast to the blowdown regime (Fig. 4) all temperatures including 
that of the LH2/LG2 interface, Tf, rise because of continuing heat absorption. Here again,the interface 
temperature is less than that of LH2 due to evaporation. The temperature stratification visualization is shown 
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Fig. 10 Tank in a blowdown regime at microgravity: 
boundary layer can be seen at t > 4000s.) 
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temperature stratification. (A slow development of the 
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Fig. 11 Tank in a storage regime at normal gravity without venting. 


in Fig. 12. Everything starts out very cold, and temperature stratification occurs naturally as heat is added to 
the tank wall. The ullage heats up while the liquid, due to its greater heat capacitance, stays relatively cool. 

Introduction of a vent valve in the upper part of the ullage volume dramatically changes the picture, 
as shown in Fig. 13. The valve is set to operate between lower and upper thresholds at 2.00 and 2.02 atm, 
correspondingly. When the upper threshold is reached, the valve opens; it closes as soon as the lower limit is 
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Fig. 12 Tank in a storage regime without venting at normal gravity: temperature stratification. 


achieved. The gas flow rate through the vent valve J v (g),e = ~ J valve where for the choked flow regime [ 22 ] 


\ Q 

J valve ~ A Pv(g)dv, valve ^ 


l{Vv + Pg) f 2 
\ Pv+Pg \7+l 


7 + 1 
7-1 


(48) 


Here A is a dimensionless relative valve position which assumes values between 0 and 1. For a given compo- 
sition of the GH2/GHe mixture, parameter 7 = 1.51, and the conditions for the choked flow regime [22] are 
satisfied. Also, in the present modeling, the valve’s cross-section area is set as S va i ve = 5 cm 2 . Each time the 
valve opens, a certain amount of the gas/vapor mixture from the upper layer is removed (Figs. 13b and 13c). 
Since that amount is small compared to the total ullage mass its removal does not substantially change the 
ullage pressure (Fig. 13a), which is uniform across the ullage volume and propagates very quickly (with the 
speed of sound). Yet the removed mass is comparable with that of the upper layer, where mass reduction, in 
accordance with the equation of state (38), leads to a temperature spike (Fig. 13f). A subsequent temperature 
wave can be observed propagating in the downward direction due to downward convection circulation in the 
bulk layers (visualized in Fig. 14), which diminishes in intensity as it travels downwards. As a result, on 
average the temperatures rise at a smaller rate compared to the non-vent case (Fig. 12), and the interface 
temperature is being stabilized (Fig. 13f). Note, that due to fast propagation of the pressure, the temperature 
variations in the lower layers follow that of pressure and are opposite in phase to those of the upper layers. 


IV. Conclusions 

In the present paper, the dynamical model of a cryogenic fuel tank has been presented that incorporates 
effects of (0 temperature stratification driven by the natural convection in both the liquid and ullage control 


20 




X lo 5 



(a) Ullage Pressure 



Time ( s) 

(d) LH2 Temperatures 



Time (s) 
(b) GH2 Mass 



(e) LH2 Bulk Temperatures 



0 1000 2000 3000 4000 5000 

Time (s) 


(c) GHe Mass 



0 1000 2000 3000 4000 5000 

Time (s) 


(f) GH2/GHe Temperatures 


Fig. 13 Tank in a storage regime at normal gravity with venting. 





n 




(a) t = 800 s (b) t = 900 s (c) t = 1000 s (d) t = 1100 s (e) t = 1200 s (f) t = 1300 s 

Fig. 14 Tank in a storage regime with venting at normal gravity: temperature stratification. 
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volumes and (ii) highly nonequilibrium condensation-evaporation phenomena at the interface between those 
volumes. Unlike standard full-scale computational fluid dynamics approach dealing with numerical solu- 
tion of equations in partial derivatives, the proposed low-dimensional scheme reduces modeling to solution 
of a set of ordinary differential equations that account for basic conservation laws applied to the properly 
chosen control volumes. Those equations are demonstrated to successfully describe complex temperature 
stratification effects for different modes of the liquid hydrogen fuel tank operation. 
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Firstly, it has been shown that the above effects are consistent with the description of the system in the 
framework of the lumped-parameter method. Secondly, time evolution of the temperature data collected 
by the fixed-elevation sensors can be described and analyzed in the framework of the developed approach. 
Thirdly, complex effects such as temperature waves caused by a vent valve operating between given pressure 
thresholds can be found and understood by means of the simplified computational approach. Fourthly, prac- 
tically important and theoretically complex effects of normal, elevated and micro-gravity can be taken into 
consideration including their role in the heat exchange between the tank walls and the cryogenic fuel con- 
tent. In particular, it has been demonstrated that the increased (decreased) gravity will substantially promote 
(suppress) such an exchange. 


Appendix A: Uniform Heat Flux 

In the above consideration, it is assumed that the tank is described by the uniform temperature of its body, 
T w . In the case of the uniform heat flux ( q w = const), the local (average) Raleigh numbers Ra x {Ran) that 
define the heat transfer coefficients (9), (22) as well as the velocity (14) in, and the thickness (15) of, the 
boundary layers (via the Grashof numbers (12)), must be substituted by a modified Raleigh number [15]. It 
can be found by expressing Ra in terms of q w by eliminating T w — T v (iy 

gp 2 c p f3q w x A 


Rat = Ra x Nu x = 




(49) 


so that the heat transfer coefficients in (9) are to be substituted with, correspondingly, the modified heat 
transfer coefficients [15]: 


K Ar K 

OLi = — NUi = — 
Xj X j 


and 


K K 

a L = j^Nuh = — 


0.631(f?.a*T-) 1/5 , 10 5 < Ra* x < 10 9 
0.241(i?a*^') 1 / 4 , 10 9 < Ra* x < 10 11 


0.726(Ra* H ^) 1/5 , 10 5 < Ra* H < 10 9 


(50) 


(51) 


0.241(i?a^T') 1 / 4 , 10 9 < Ra* H < 10 11 

where T is given by (10). Given the uniform heat flux q w the local temperature at the wall surface can be 
found from 


T =T 


v(l) 


Qw 


q w x 


a V< ' 1 ' > Nu x k 


(52) 
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where the Nu x number can be found from (50). 


Appendix B: Lumped-parameter Approximation 

Here is described the lumped-parameter approximation for each of the three CVs: liquid, vapor and film, 
similar to that presented in [1] by ascribing a common temperature to each of those CVs. In that case the 
mass conservation equation for, e.g., the liquid CV can be obtained by summing up (26), (28), (31), and (35) 
over all liquid elements thus yielding 

tiq — Jlv (53) 

Here the interface flow rate can be found by solving (7.1) together with the total-mass conservation for 
the vapor CV 


TTl v — Jlv “I” Jv,e 

(54) 

rilg = Jg,e 

(55) 


and the energy conservation for the liquid, vapor and film CVs (see also (47)), all of which can be obtained 
in a similar way,: 

Qi,e + Qi + W + Jivhis + Ji , e (hi + vf e /2) = —miui (56) 

’ at 

Qv,e Qv if Jlvhys ~f Jv,e(hvV v , e /2) Jg e {hg 4- Vg e j2ij — -^(TTl v U v -f THgUg') (57) 
complemented with the expression relating the liquid and vapor volumes 

V = V v + V l = V v + m l /p l (58) 

and the film (39) and vapor(gas) equations of state 

P y {g ) Pv(g) ^v(g) r Pv ( 59 ) 

Given the external mass flow rates, J v>e , Ji >e , Eqs (53) - (59), (39), and (47) allow one to find the rate of 
the total mass change rhi (or Vi) for the liquid CV, the rate of the volume change V v of the vapor CV, as 
well as the “average” temperatures T/, T v as functions of time. In this step, the local convection heat transfer 
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coefficients for the vertical surfaces should be substituted by their average values [15] 


0.68 + 0.67 (Ra H ^) 1/A , 10 5 < Ra H < 10 9 

a H = 77 — < (60) 

M l(v) 

I 0.15(Ra H ^) 1/3 , 10 9 < Ran < 10 11 

where Hu v ) is the height of the liquid/vapor in contact with a vertical wall. 


Appendix C: Model Parameters 

Model parameter values are given in Table 1. Here, d w refers to the tank wall thickness. 

Table 1 Parameter Values 

Component Parameter Values 

Hydrogen and Helium T c = 33.2 K, p c = 1.315 x 10 6 Pa, A = 5, p L = 71.1 kg/m 3 , c L = 9450 J/kg/K, 
k l = 0.0984 W/m/K, p = 3.4 x 10~ 6 Pa s, R v = 4124 J/kg/K, 
c v , v = 6490 J/kg/K, 7 = 5/3, k v = 0.0166 W/m/K, R g = 2077 J/kg/K. 

K g = 0.0262 W/m/K, c v , g = 3121 J/kg/K, c P , 9 = 5193 J/kg/K 
Tank n = k = 10. R = 4.21 m, H = 29.56 nr, c w = 1 x 5 x 10 2 J/K, p w = 2700 kg/m 3 , d w = 0.1 m, 

*5 valve — 0.05 cm" 
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